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Abstract 

In the spin-boson model, a continued fraction form is proposed to systematically resum high-order 
quantum kinetic expansion (QKE) rate kernels, accounting for the bath relaxation effect beyond 
the second-order perturbation. In particular, the analytical expression of the sixth-order QKE rate 
kernel is derived for resummation. With higher-order correction terms systematically extracted 
from higher-order rate kernels, the resummed quantum kinetic expansion (RQKE) approach in the 
continued fraction form extends the Pade approximation and can fully recover the exact quantum 
dynamics as the expansion order increases. 
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I. INTRODUCTION 


In a quantum dynamic process, the interaction between the system and bath leads to 
irreversible energy relaxation and decoherence of the quantum system. The resulting quan- 


Q. 


turn dissipation can lead to rich quantum phenomena, e.g., quantum phase transition 
The spin-boson (Caldeira-Leggett) model is a simple but fundamental quantum system, 
which can be used to interpret the quantum tunneling and localization in macroscopic sys¬ 
tems 2, 3]- Gate operations in quantum computation and quantum information are sim¬ 
ulated by quantum dissipative dynamics of multiple spin-boson models, where each qubit 
is equivalent to an individual spin { 4 ]. In the study of quantum transport, a fundamental 
question is to understand the transport process from a donor to an acceptor in the two- 


site system [5|J. In the lowest order, the transfer rate is estimated using Fermi’s golden 
rule (FGR), proportional to the square of the site-site coupling strength. This second-order 
transfer rate is expressed as the Forster theory in energy transfer |6J and as the Marcus 
theory in electron transfer [7]. The non-Markovian relaxation of the surrounding bath can 
significantly slow down the transfer process compared to the second-order prediction 8-17]. 


On the other hand, the transfer rate can be optimized at an intermediate dissipation strength 



As a simple quantum model, the spin-boson model (or the equivalent two-site system) 
is a benchmark system for the study of quantum dynamic methodologies. In addition to 
the sophisticated Feynman-Vernon influence functional 26], a straightforward app roach of 


quantum dissipation is to apply the Nakajima-Zwanzig projection operator [27], 28]. In 
the lowest second order, we obtain various approximate dynamic equations from different 
perturbed terms, e.g., the Redheld equation from the system-bath interaction 29], and 
the FGR rate from the site-site coupling. The noninteracting-blip approximation (NIBA) 
extends the FGR rate to a time-nonlocal description of the detailed time evolution js|. 
To improve the NIBA prediction, the variational polaron method is a modified second- 
order perturbation where the perturbed term is self-consistently determined from equilibrium 
distribution |30|, [31]. The variational polaron method is more reliable in the unbiased two- 
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site system with a relatively fast bath. A more systematic approach beyond the second-order 
perturbation is to include higher-order corrections of perturbed terms, as in the quantum 
kinetic expansion (QKE) approach 12417]. In our recent paper 16], the higher-order QKE 
of the site-site coupling is obtained using an indirect projection operator technique for a 
general multi-site system. In the two-site system, all the higher-order QKE corrections arise 
from the bath relaxation effect, whereas in the multi-site system, the higher-order QKE 
corrections also include quantum interference effects. 

A key theoretical concern in the QKE approach is the resummation technique of higher- 
order rate kernels, due to two essential reasons. The analytical and numerical difficulties 
quickly increase as the expansion order increases. More crucially, the QKE rate kernels can 
converge slowly and become divergent as the site-site coupling increases. An appropriate 
resummation technique can partially include corrections of all the orders using one or a few 
higher-order QKE rate kernels, and avoid the divergence of large site-site couplings. For 

echniques are the Pade approxima- 


the lowest-order correction, two typical resummation 

tion [l2, 13] and the Landau-Zener approximation ]32]. With a factor of 2 difference, the 


Landau-Zener approximation is not reliable in the strong dissipation limit, compared to the 


Pade approximation. In a recent paper 


33], a modified resummation approach is proposed 


with an optimization according to the equilibrium distribution. However, any resummation 
techniques in the lowest order cannot fully account for the extra knowledge of higher-order 
QKE rate kernels, and its prediction deviates significantly from the exact quantum dynamics 
at some point. 

Therefore, a more general resummation technique is required to systematically include 
corrections from higher-order rate kernels. In Ref. |34j, a generalized Pade approximation is 
developed, which is complicated in its mathematical formulation and practical application. 
Instead, we will extend the physical factorization scheme in the Pade approximation to the 
higher-order QKE rate kernels and obtain a simple continued fraction form, which leads 
to a systematic resummed quantum kinetic expansion (RQKE) method. In Section [TT1 the 
derivation of the QKE approach in the two-site system (the spin-boson model) is briefly 
reviewed. The time-integrated QKE rates of the first three orders are numerically computed 
in a quantum Debye bath. In Section IIII1 the continued fraction resummation form is 
developed, and the RQKE rates are numerically compared with the exact results of both 
unbiased and biased systems. In this paper, all the exact quantities are obtained using the 
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hierarchy equation method 354381]. In Section llVl the RQKE rate kernels are used to predict 
the detailed population evolution, and are calibrated with the exact result. In Section El the 


temperature-dependent equilibrium population is calculated using t 


is also compared with the exact stochastic path integral result 
summarize our studies. 


39 


re RQKE rates, which 


40]. In Section EH we 


II. QUANTUM KINETIC EXPANSION IN A TWO-SITE SYSTEM 


In this section, we briefly review the quantum kinetic expansion (QKE) approach in 
Ref. [h]]. With respect to the single excitation manifold, the bare Hamiltonian of a multi- 
site system is given by H$ = J2 n £ n\ n ){ n \ + Ylm^m Jmn\m)(n\, where | n) represents a single¬ 
excitation quantum state localized at site n, e n is the excitation energy of site n, and J mn 
is the site-site coupling strength between sites m and n. The bare Hamiltonian of the 
surrounding environment is given by Hb■ The system-bath interaction Hsb is considered to 
be localized at each site n, H SB = HsB-nln)(n\. In the site basis representation {|n)}, 
the total Hamiltonian is written as 


tftot ^ ^ H n 177.) (Tl | T ^ ' Jmn |m) (jl |, (1) 

n nm(n^m) 

with H n = e n + Hb + Hsb;ti- Here the simplest two-site system coupled with a harmonic bath 
can be mapped to the standard spin-boson model. The time evolution of the total density 
matrix p tot (t) follows the Liouville equation, d t p tot (t) = -iC tot p tot (t), with C tot = [H tot , • • •]• 
Throughout this paper, the reduced Planck constant h is treated as a unit. Following the 
separation of population and coherence components, the total Liouville superoperator is 
formally expressed as a block matrix, 


( £p £pc \ 

y £cp J 


( 2 ) 


where the subscripts P and C denote system population and coherence, respectively. In the 
two-site system, the diagonal part of C to t is fully dependent on the diagonal Hamiltonian 
elements H n , while the off-diagonal part of £ to t arises from the site-site coupling J . Sub¬ 
sequently, we define the partial time propagation superoperators, Up(t) = exp (—i£ P f) and 
U c (t) = exp(— iCct), which can be interpreted as Green’s functions in the Liouville space. 
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An indirect projection operator approach is applied in Ref. 16] to derive the closed time 
evolution equation of the reduced system population P(t). The initial condition is required 
to be a local equilibrium state, ptot(O) = YJ n PnP*Bn\ n )( n \i where p n is the initial population 
of site n, and pg q n oc exp (—/3H n ) is the local Boltzman density of bath. The final time 
evolution equation of P{t ) follows a time-nonlocal convolution form, 

P{t) — — f drlC(t — r)P(r). (3) 

Jo 

The rate kernel /C(f) is derived as an expansion form of the site-site coupling J, given by 
/C = /C® + /C (3) + • • •. In the two-site system, all the odd-order terms vanish, and only the 
even-order terms remain. Here we introduce a local equilibrium population state matrix, 


p(°) = 
re q 


P B ;l 0 


0 {% 


(4) 


and its projection matrix, viq = piq }Ttb{, where Tr B {- • • } is the partial trace over bath 
degrees of freedom. The 2/c-th QKE rate kernel is explicitly given by 

£ (2fc) (>2,T3, • • • ,r 2fc ) 

= — (—l) fc Tr B {[7?.(r2fc)(jWp(T2fe- 1 )] 

[ft(T 2 fc_2)<5Wp(T2fc- 3 )] • (5) 

where 5Up(t) = Up(t) — Vcq is the pure dissipative propagation superoperator, vanishing in 
Markovian dynamics, and Pit) = CpcUcJt)C cp is the population-to-population transition 
superoperator. Thus, high-order {k > 2) QKE rate kernels reflect dynamics of population 
fluctuation around the local equilibrium state due to the bath relaxation effect of 5lip{t). 
Equation (J^D is equivalent to the previous expression of Eq. (15) in Ref. 16], but in a more 
concise form. The Feynman diagram technique is applied to visualize these quantum rate 
kernels in Fig. |TJ which is also simplified in notation compared to previous diagrams in 
Ref. 1|J. In detail, each initial and final numbered circle represents a local equilibrium 
population state, p‘p. ri \'n)(n\, at the corresponding site n(= 1,2). Each intermediate dashed 
circle represents the dissipative propagation [§Up(t)\ n of a system-bath entangled population 


state, pp-nif) — [ptot(j)]nn- Unlike the notation in Ref. 16], each arrowed line represents a 


population-to-population transition [P(t)\ mni as a density flow from population to coherence 
and back to population, [JZ(t)] mn = \J m n\ 2 \JJc-,mn{t) + Uc-nm(t)\- 
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FIG. 1: The Feynman diagrams of the second- (a), fourth- (b), and sixth-order (c) quantum rate 
kernels in the two-site system (the spin-boson model). The explicit interpretation of each symbol 
is provided in text. 

The formal expression of /C^ 2 fc ^(r 2 , 13 , • • • , T 2 k) in Eq. (J5]) is derived for an arbitrary en¬ 
vironment, beyond the spin-boson model. Next we assume that the bath is harmonic and 
Hsb follows a bilinear form. With the creation (a) 1 ") and annihilation (a*) operators for the 
it\i harmonic oscillator, the bath-coupled Hamiltonian at local site n reads 

H n = £ n + ^2 ^4 a i + ^2 U i x ni {<4 + Osi) , ( 6 ) 

i i 

where the coefficient x ni denotes the system-bath coupling strength reduced by the frequency 
Ui of the ith harmonic oscillator. The QKE rate kernels in Eq. ([5]) are transformed into the 
time correlation functions of the displacement operator, G n = exp [ Yhi x m{ ( 4 ~ a;)], which 
can be obtained by the cumulant expansion. If the bath coupling is identical at each system 
site, the explicit expression of the second-order rate kernel reads 

K'rnn(j^m)44 2|t/ mn | Re exp{ [Kmnh T S mn Q^T 2 )]}, (7) 

where £ nm = e n — £ m is the modified site excitation energy detuning with i n = £ n — 
'52i U} i x nn an d ^ ie coefficient s mn arises from the site-site ‘spatial’ correlation. For the 
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(a) 


(b) 



FIG. 2: The normalized time-integrated forward transfer rate expansions of the first three orders 

(a) in the unbiased system with £12 = 0, and (b) in the biased system with £12 = 100 cm -1 . 

The Deybe frequency of the coupled bath is = 100 fs, and the temperature is T = 300 K. 

The normalization is realized by (a) k^^ D {\)/k^} D {\ = lcm” 1 ), and (b) k^^} D (X)/k^^} D . max . 

( 2 ) 

In each figure, the dotted black line is the second-order result k\^_ D , the dashed blue line is the 
fourth-order result k^_ D , and the solid red lines is the sixth-order result k^_ D . 


standard spin-boson model, a perfectly negative correlation leads to = 4, while for 

the regular energy transfer system, a 5-spatial correlation leads to = 2. Thus, the 

two-site system under the 5-spatial correlation is equivalent to the spin-boson model with a 
doubled dissipation strength (reorganization energy). The time correlation function of the 
displacement operator excluding the spatial dependence is 


g(t) = / dw[J(a;)/(n 2 ][(l — coswf) coth(/3cu/2) 

Jo 

+ i sin ut ], 


( 8 ) 


where J{oj) = (u> — u>i) is the bath spectral density. Equation (JTJ) is the same as 

the rate kernel in the NIBA approach [ 3 }], and its time integration recovers the FGR rate. In 
Ref. [uij], the fourth-order QKE rate kernel is derived for a general multi-site system. The 
simplified expression of 73, 74) for the two-site system with the 5-spatial correlation 

is provided in Appendix 0 Furthermore, we extend to the sixth-order QKE rate kernel, and 
the explicit expression of 16 terms is also shown in Appendix lAl 

Before investigating the resummation technique in next section, we numerically calculate 
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the quantum rate kernels of the first three orders. Both unbiased and biased two-site systems 
are considered with Syx = 0 and 100 cm -1 . To be compared with the calculation of the 
hierarchy equation I35l-l38|. a quantum bath with the Debye spectral density is applied, 
given by 




2 ’ 


(9) 


\ 7T J U 2 + U) D 

where @(ca) is the Heaviside step function of oj, A is the reorganization energy, and ojjj is the 
Debye frequency. For simplicity, we introduce the high-temperature approximation, leading 
to 


9(t) 


2A 


\t\- 


1 _ e ~ UJ D |£| 

OJD 


+ ?Sign(t)A 


X _ e -u D \t\ 
U}D 


( 10 ) 


where Sign(f) is the sign function of t. In our calculation, the Debye frequency is = 100 
fs, and the temperature is T = 300 K. We focus on the time-integration of rate kernels, 
^( 2 fc) _ J o °° dr 2 • • •/ 0 °° dr 2 A,-/C^ 2fc ^(r 2 , • • • , r 2 fc), which can be viewed as the time-integrated 
effective rate matrix, especially for over-damped dynamics. Since the 2/c-t.h rate kernel is 
proportional to the 2/c-th power of the site-site coupling J, we normalize effective rates 
to remove the J-dependence. The normalization is over the maximum value /Cmal for the 
biased system, and over the value of the minimum reorganization energy (A = 1 cm -1 ) for 
the unbiased system. Due to the heavy computational duty in a multi-time integration, the 
Monte Carlo simulation of 10 12 samples is applied to the calculation of /C^ for convergence. 
Figure [2] presents the numerical results of the forward transfer rate expansions from 

the donor site 1 to the acceptor site 2, which will be used for the resummation technique 
in next section. We find that k^_ D monotonically decreases with the reorganization energy 
A in the unbiased system, whereas kr^_ D is maximized in an intermediate value of A in the 
biased system. 


III. RESUMMATION OF QKE RATE KERNELS IN A CONTINUED FRACTION 
FORM 

In the previous section, we present the explicit expansion forms of rate kernels in the two- 
site system (the spin-boson model) using the QKE approach. For a small site-site coupling 
strength, the full quantum kinetic rate kernel can be obtained straightforwardly as the sum 












FIG. 3: The Feynman diagrams of the fourth- (a) and sixth-order (b) quantum rate kernels in the 
two-site system (the spin-boson model) under the Pade approximation. The matrix factorization 
is realized by inserting vertical dashed lines. The other symbols are the same as those in Fig. |T] 


of /C (2fc ) up to a converged expansion order. For a large site-site coupling strength, this simple 
summation cannot be applied since /C^ 2fc) diverges as the expansion order increases. Instead, 
a resummation technique is required for a converged result, with one or more high-order 
corrections of (k > 2). For the leading-order QKE correction various resummation 


methods, e.g., the Pade approximation 


12] and the Landau-Zener approximation [32], have 


been well discussed previously. Although these approximations can significantly improve the 


second-order prediction of the NIB A approach jl2-ll7j], a systematic resummation approach 


is still required to include higher-order corrections and recover the exact quantum dynamics. 

We revisit the Pade approximation in Ref. 12] to show its physical interpretation, which 
will used for a generalized resummation technique. As mentioned in previous section, the 
pure dissipation of population, 5Up(t), vanishes in Markovian dynamics. For a fast relaxing 
bath with a weak non-Markovian feature, or alternatively in the strong dissipation regime 
where the system transport is slow but Markovian, an approximate time separation can be 
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expected in the high-order QKE rate kernels. For the leading-order correction /Q 4 ) ( 72 , 73 , 74 ), 
this approximation is realized mathematically by inserting a reduced population projection 
Vp before the action of 5Up{r 3 ) [12| . In the reduced population subspace, Vp is explicitly 
written as 


V P 


( PB q }Tr B { 0 \ 

{ 0 p^}Tr B { ) 


( 11 ) 


where pg q oc exp(— (3H B ) is the bare bath equilibrium distribution. Equation (fTTjl results 
in two identities, VpVeq = Vp and VeqVp = Veq ■ As a result, the fourth-order QKE rate 
kernel is factorized into 


£ (4) ( r 2, 73,74) « 5 (2) (r 3 ,r 4 )/C (2) (r2), (12) 

with ( 73 , 74 ) = —Tt b {TZ(t4 : )5Up(t 3 )Vp}. The matrix factorization can be applied to all 
the higher-order corrections, giving 

/C (2A) (r 2 , • ■ ■T 2 fc) 

~ S (2) (t2/c-i, r 2 k) ■ ■ ■ S (2) (r 3 , r 4 )/C (2) (r 2 ). ( 13 ) 


Figure [3] presents the Feynman diagrams of /Q 4 )(t 2 , r 3 , r 4 ) and /C < - 6 ' ) (r 2 , • ■ -tq) after the matrix 


factorization. With the introduction of the 
correction term S( 2 )(;Q of lC^\z) becomes 

r( 4 ) ( z ) — X 

'Wesuml^v - L ' 


japlace ^-transform, the resummation using the 


12 ] 


7 ( 2 ), 


1 -1 


M K <2| M. 


(14) 


where X is an identity matrix. By expanding Eq. (fT4l) in the 2x2 matrix form, we recover the 
regular Pade approximation for both forward (&resum-A-<-.D(' 2 0) an< l backward (^resum-rn-A^)) 
transfer rate kernels. 

Next we can extend to higher-order corrections with a generalized factorization technique. 
Following the definition of S( 2 )(t 3 ,t 4 ) to higher-orders, we introduce another expansion 
series, 


X (2fc) (r 3 , • • • ,t 2 *) 

= (—l) fc Tr B {[7^(r2fc)5Wp(r 2 fc-i)] ■ ■ ■ [K(u)6Up(t 3 )]V p }, 


(15) 


which is essential for the QKE in the system-bath separated initial condition 4l|. For 
the sixth-order QKE rate kernel, a more accurate matrix factorization is changed to 
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(a) (b) 



FIG. 4: The RQKE forward transfer rate from the continued fraction form compared with its exact 
value from the hierarchy equation in the unbiased system with (a) J = 40 cm” 1 , and (b) J = 100 
cm” 1 . In each figure, the black dotted line is the second-order Forster rate, the blue dashed line 
is the lowest fourth-order RQKE rate (the Pade approximation), the red solid line is the next 
sixth-order RQKE rate, and the green circles are the exact result. The parameters of u>d and T 
are the same as in Fig. [2j 


/C (6 )(t 2 ,''' , tq) ~ S( 4 )(t 3 ,--- , T 6 )/C (2 )(t 2 ). Similar to the cumulant expansion, the ‘real’ 
fourth-order correction needs to exclude the contribution of S®, 

5S (4) (t 3 ,--- ,t 6 ) 

= H (4) (t 3 , • • • ,r 6 ) - E (2 ) (t5,t 6 )H (2) (t3,t 4 ). (16) 

All the other higher-order QKE rate kernels are subsequently factorized using Si 2 ) and 5EA 4 '. 
For conciseness, we introduce the difference of 5EA 4 ' relative to Si 2 ), which is defined in the 
Laplace z-space as 

E< 4 >( 2 )=A 4 ( z )S< 2 >(z). (17) 


Here the expansion index 4 is assigned as a subscript since A 4 (z) is in the same J-expansion 
order as A 2 (z) = 5^ 2 ^(z). The approximate full quantum rate kernel resummed from A 2 (z) 
and A 4 (z) is derived in a continued fraction form, 


£( 6 ) 


(*) 



A 4(2) 



A 2(2) 


/C (2) (^). 


(18) 
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The above factorization scheme can be straightforwardly to an arbitrary expansion order, 
which defines the general correction term, §’E^‘ 2k \z) = A. 2 k- 2 (z ) ■ ■ • A 2 (z), and gives rise to 
the general matrix continued fraction form. 

The separation of higher-order QKE rate kernels depicted in Fig. [3] requires modifica¬ 
tions when the non-Markovian dynamics is not weak. The dynamic coupling between 
H( 2fe — 2 )(r 3 , • • ■ , T 2 fc) and KL^ 2 \t 2 ) needs to be include, beyond the matrix factorization, 
/C (2fc) (T 2 , • • • , T 2 fc) ~ 5 (2fc ~ 2) (r 3 , • • • , T 2 fc)/C (2) (r 2 ). However, this difficulty can be circum¬ 
vented using the scalar continued fraction form for each element of the rate kernel. Mathe¬ 
matically, a regular function can be re-expressed in the continued fraction form, by matching 
its Taylor expansion series. Thus, we propose the scalar continued fraction resummation 
form for the forward rate kernel, 

£m -A^d(z) = -7- ~k { A ( 19 ) 

$2 ;A<-d(z) 

1 + - 

1 + ^2fc-2;A<-D(^) 

where the correction terms are matching the QKE forward rate kernels k^_ D {z) term by 
term, given by 

S 2 ; a<-d(z) = -k%l D (z)/k ( Jl D (z), (20) 

S 4 ;A^d(z) = -S 2 ;A^d(z) - kf^ D (z)/kf^ D (z) , (21) 


The same approach is applied to resum the backward rate kernel Equa¬ 

tions (USD-flnD provides the basic construction of the resummed quantum kinetic expansion 
(RQKE) method. To be consistent, the expansion order of the RQKE is defined by the 
power of the site-site coupling strength in the highest-order QKE rate kernel considered. 
Compared to the generalized Pade approximation in Ref. 34], the continued fraction can 
also be expanded into a rational polynomial form, while the correction terms in the RQKE 
method are more straightforwardly obtained without an additional basis expansion. In addi¬ 
tion, as the resummation order 2k increases, all the lower-order correction terms S 2 j(<k-i)(z) 
are not affected, which makes the continued fraction form a more systematic approach. 

To verify the reliability of the continued fraction form, we use the results of the first three 
order effective rate expansions in Section [TT] to obtain the RQKE rates fcSm = krcsum(z = 0), 
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(a) (b) 



FIG. 5: The RQKE forward transfer rate from the continued fraction form compared with its exact 
value from the hierarchy equation in the biased system, £12 = 100 cm -1 , with (a) J = 20 cm -1 , 
and (b) J = 100 cm -1 . In each figure, the black dotted line is the second-order Forster rate, the 
blue dashed line is the lowest fourth-order RQKE rate (the Pade approximation), the red solid line 
is the next sixth-order RQKE rate, and the green circles are the exact result. The parameters of 
oj£) and T are the same as in Fig. [2j 


which are compared with the exact full quantum rates k ex ac t from the hierarchy equation. 
In Ref. jli|, fcexact is calculated under a system-bath separated initial condition, different 
from the presumption of the local equilibrium population state in the QKE approach. The 
accurate value of k exact is re-calculated, following the rigorous expression in Ref. [4l|. With 
the same equilibrium population, the results of k exact under these two initial conditions are 
proportional to each other 41]. The results of k^, kitlum , fcresum, and k ex ac t for the forward 
transport process form the donor site 1 to the acceptor site 2 are plotted in Figs. [HandO 
For the unbiased system (ei 2 = 0), two site-site coupling strengths, J = 40 and 100 cm -1 
are considered; for the biased system (ei 2 = 100 cm -1 ), two site-site coupling strengths, 
J = 20 and 100 cm -1 are considered. For the two small site-site coupling strengths, J = 40 
cm -1 and e 12 = 0 in Fig. Uk, and J = 20 cm -1 and e 12 = 100 cm -1 in Fig. [5^, the QKE rate 
kernels converge with the expansion order. The lowest fourth-order RQKE rate k^l um . A ^ D 
improves the second-order FGR rate and predict k exaC f } A<-D accurately in the whole range 
of the reorganization energies, 1 cm -1 < A < 1000 cm -1 . For the large coupling strength of 
J = 100 cm -1 in Figs. SJo and[5]o, the QKE rate kernels diverge with the expansion order. In 
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the unbiased system, kltlum-A^o improves k^_ D mainly in the large-A regime. In the biased 
system, k^l um . A ^ D largely improves except for a small deviation in the intermediate- 

A regime. In Fig. 5b of Ref. 16], the difference between k^ sum . A< ^ D and /c eX act ; n<-D in the 
large-A regime is due to an inconsistent initial condition in the hierarchy equation. For 
both cases, the next sixth-order RQKE rate k^l um . A ^ D agrees perfectly with k exa , ct -A<^D in 
the whole A regime. Our numerical calculations demonstrate that the RQKE rate from the 
continued fraction form can systematically approach to the exact value, and the number of 
necessary correction terms gradually increase with the site-site coupling strength. 


IV. TIME-CONVOLUTED QUANTUM KINETICS 


The continued fraction form of the bath relaxation effect is verified by the convergence 
of the resummed effective rate toward the exact value. In this section, we will further 
demonstrate the accuracy of the continued fraction in predicting the detailed time evolution 
of site population. 

All the high-order QKE rate kernels can be derived explicitly, using the cumulant expan¬ 
sion for the multi-time correlation function of the displacement operator. The time evolution 
of reduced site population P(t ) is subsequently solved by the convoluted equation in the time 
f-space, or equivalently by the matrix inversion in the Laplace ;;-space. The computational 
cost of both methods is often very high. Instead, we re-express the QKE rate kernels in a 


matrix 
in Ref. 


ormalism 41]. The general 2/c-th QKE rate kernel in the Laplace £-space is derived 


4l| as 


K [2k \z) = ~T P [n(z)SU i0 \z)} k - 1 iZ{z)V^, (22) 

with iZ(z) = W (1 ^W ( °\^)W (1) . Here each matrix is defined in an expanded basis set of 
relevant dynamic variables and can be mapped to a superoperator in Section Q4] Specif¬ 
ically, the mapping of two projection matrices are Vp Vp and Vcq ^ V^. The two 
interaction Liouville superoperators are combined together and mapped to a perturbed tran¬ 
sition rate matrix, {i£pc,i£cp} O W (1) . The two unperturbed time propagation superop- 
eraotrs are also combined together and mapped to an unperturbed pure dissipative matrix, 
{5Up(z)Mc(z)} 8U [0 \z). 

For over-damped quantum dynamics in the two-site system, the time evolution of site pop- 
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FIG. 6: The time evolution of the donor population in the unbiased system with £12 = 0. The 
same quantum Deybe bath in Fig. [2] is applied. The conditions of the four figures are (a) J = 20 
cm” 1 and A = 4 cm” 1 , (b) J = 20 cm” 1 and A = 12 cm” 1 , (c) J = 100 cm” 1 and A = 4 cm” 1 , and 
(d) J = 100 cm” 1 and A = 12 cm” 1 . In each figure, the dashed line is from the exact dynamics, 
and the dashed-dotted line is from the lowest-order resumed kernel KAttL m (z) . In the top two 
figures, the dotted lines are from the second-order kernel JC^(z). In the bottom two figures, the 
solid lines from higher-order resumed rate kernels fully recover the results of the exact dynamics 
using Eq. (|23l) and coincide with the dashed lines. 

illation is close to a single exponential decaying function (Markovian behavior), which can 
be described by the time-integrated effective rate. To illustrate the relevant non-Markovian 
behavior, we focus on small and intermediate reorganization energies with under-damped 
dynamics. In our two-site system, we choose two typical reorganization energies, A = 4 
and 12 cm” 1 , for each system condition (£12 and J) in Figs. [4] and El The exact time 
evolution of site population, P exa , c t-,i{.t), is solved by the hierarchy equation using the local 
equilibrium population state initially at the donor site 1. Next we re-calculate the site pop¬ 
ulation P exac t-i(z) in the Laplace 2 -space, and obtain a new estimation of the time evolution, 
Pg xac t i(^) = LT” 1 [P eX act ; i(^)], using the inverse Laplace transform, LT _1 [- ■ •]. The two time 
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FIG. 7: The time evolution of the donor population in the biased system with S 12 = 100 cm -1 . The 
same quantum Deybe bath in Fig. [2] is applied. The conditions of the four figures are (a) J = 20 
cm -1 and A = 4 cm” 1 , (b) J = 20 cm” 1 and A = 12 cm” 1 , (c) J = 100 cm” 1 and A = 4 cm” 1 , and 
(d) J = 100 cm” 1 and A = 12 cm” 1 . In each figure, the dashed line is from the exact dynamics, 
and the dashed-dotted line is from the lowest-order resumed kernel /citL ™ (z) . In the top two 
figures, the dotted lines are from the second-order kernel )C^ 2 \z). In the bottom two figures, the 
solid lines from higher-order RQKE rate kernels fully recover the results of the exact dynamics 
using Eq. (l23l) and coincide with the dashed lines. 


evolution predictions, -P e xact ; i(£) and Pexact-iW are found to be identical, confirming the re¬ 
liability of the numerical inverse Laplace transform. In onr model system, Eq. (12 2 p is also 
numerically solved using the hierarchy equation approach 4lJ. The estimation of the site 
population from the 2/c-th order RQKE rate kernel is written as 


p£LM = LT 


-1 


k 


(2k) 

resum;A<— D 


» 


Z [ Z + ffim;A^ d( Z ) + KtsL^A^)}. 


(2k) 


(23) 


We apply the same two-site system with the same quantum Debye bath with c= 100 
fs and T = 300 K in previous two sections. The comparison between P^\t), P^j nr i (t) and 
-Pexact ; i(fO is organized in Figs. [HI and [TJ where P, (t) is the second-order NIBA prediction. 
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In the unbiased system (e = 0) with the small site-site coupling (J = 40 cm -1 ), P± 2 \t) 
is close to the exact time evolution P e xact;i(£) with a small deviation. The lowest fourth- 
order RQKE rate kernel, X)resum(^), further improves P,- 2 " 1 (t) and provides almost identical 
results of -P e xact ; i(Q for the two reorganization energies. As the site-site coupling is increased 
to J = 100 cm -1 , Prctum-iW improved from the NIB A prediction also deviates from the 
exact result P e xact ; i(^)- We find that Prcsum iW gradually approaches to P e xact ; i(^) as the 
J-expansion order 2k increases in the continued fraction form. As shown in Figs. [6t and d, 
■^resum-i (0 an d -^resum-i(^) f rom the eighth- and sixth-order RQKE rate kernels fully recover 
Pexact; i (Q for A = 4 and 12 cm -1 , respectively. In the biased system (e = 100 cm -1 ) with 
the small site-site coupling (J = 20 cm -1 ), P[~\t) clearly deviates from P exa ct;i (t), while 
Piesum- 1(0 from the Pade approximation becomes almost identical to P e xact ; i(Q for the two 
values of A. Although the time-integrated rate fcresum is very close to the exact value k exact in 
the small- A regime, the prediction of -P r p 4 j 1TTV 1 (t) is no longer reliable for the strong site-site 
coupling (J = 100 cm -1 ). Similarly, we extend the continued fraction form to higher orders, 
and Presum-iW fr° m the tenth-order RQKE rate kernel fully recovers P exact ;i (t) for A = 4 
and 12 cm -1 . Thus, the exact quantum dynamics can be fully predicted by the RQKE rate 
kernels in the continued fraction form. The convergence order of the continued fraction for 
the detailed time evolution in general increases as the reorganization energy decreases. Since 
the equilibrium population in the unbiased system is unchanged with the system and bath 
parameters, the convergence order is usually smaller than that in the biased system. 

V. TEMPERATURE DEPENDENCE OF THE QUANTUM EQUILIBRIUM POP¬ 
ULATION 

In this section, we will further demonstrate the accuracy of the continued fraction in 
predicting the temperature dependence of quantum equilibrium population. 

In the original matrix continued fraction form, the expansion (z) from the factoriza¬ 
tion scheme on the high-order QKE rate kernels leads to the same correction terms for both 
forward and backward transfer rate kernels, i.e., 5 2 k-,A<^D{z) — $ 2 k;D<-A(z)■ The ratio of the 
two time-integrated RQKE rates, k^l m . A ^ D /k^l nvD ^ A , is unchanged as the resummation 
order increases. The equilibrium population is always the same as the classical Boltzmann 
distribution of the FGR prediction, P eq -n oc exp(— j3e n ), which is only valid at high tempera- 
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FIG. 8: The equilibrium donor population versus the temperature. The solid line with up-triangle 
symbols is the result of Fermi’s golden rule rate. The solid line with diamond symbols is the result 
of the lowest fourth-order RQKE rates using the Pade approximation. The circle symbols represent 
the result of the next sixth-order RQKE rates. The solid line without symbols is the exact result 
from the hierarchy equation. The parameters are £12 = 100 cm -1 , J = 100 cm -1 , A = 100 cm -1 , 
and = 100 fs. 


tures. In our modified scalar continued fraction form, the correction terms of the forward and 
backward rate kernels are determined independently, which allows 7 ^ ^2 

Consequently, the equilibrium population predicted by the RQKE rate can deviate from the 
classical Boltzmann distribution an d appro ach to the exact quantum Boltzmann distribu¬ 
tion, P cq . n oc [Tr B {exp(-/3# tot )}] nn 33,140]. 

As a verification, we extend our previous study at a high temperature T = 300 K to lower 
temperatures. Since the equilibrium population is always one half in the unbiased system, 
we only consider the biased system, £ 12 = 100 cm -1 with J = 100 cm -1 and A = 100 cm -1 . 

(2k) 

The 2/c-th order prediction of the donor equilibrium population P eq;r g Sum;1 is obtained using 
the time-integrated RQKE rates, 


p(2fc) 

■^eq; resum; 1 


mesum;Z) <—A 


k 


( 2*0 


+ h 


(2k) 


(24) 


u:esum;A-<— D ' ^resum ;!)•<—A 

The full expression of the time correlation function g(t) is applied in the calculation of the 
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QKE rate kernels, without the high-temperature approximation. Similarly, 10 12 Monte Carlo 


samples are simulated for an accurate estimation of k 


( 6 ) 


The hierarchy equation with 


the Matsubara frequency summation is used to obtain the exact equilibrium population, 


which is numerically the same as the result of the stochastic path integral 


|39 


40j . Our 


numerical calculation shows that each correction term S‘ 2 j(<k-i) is different for the forward 
and backward rates, and the deviation increases as temperature decreases. As shown in 


Fig. El the RQKE rates systematically improves the prediction of P 


(2k) 

eq;resum;l 


from the second- 


order FGR result to the exact result. With specihc parameters in our calculation, the 
sixth-order RQKE rates provide an excellent prediction of the exact result over the whole 
temperature range (100 K < T < 300 K). With more correction terms included, we expect 
that the scalar continued fraction resummation can be straightforwardly extended to lower 
temperatures. 


VI. SUMMARY 

In this paper, we extend our previous study of the quantum kinetic expansion (QKE) 
approach in the two-site system (the spin-boson model). The factorization scheme for the 
high-order QKE rate kernels in the weak non-Markovian dynamics leads to the matrix 
continued fraction form for the resummation technique of QKE rate kernels. To be valid in 
an arbitrary condition, we further introduce the scalar continued fraction form for forward 
and backward rate kernels separately, where the correction terms are obtained by matching 
the higher-order QKE rate kernels. Consequently, a systematic resummed quantum kinetic 
expansion (RQKE) method is constructed, and the expansion order of the RQKE method 
is consistent with the highest order of the QKE rate kernel. To the lowest fourth-order, 
the continued fraction form recovers the Pade approximation, while the higher-order RQKE 
correction terms represent the additional bath relaxation effects. As shown by numerical 
calculations in this paper, the prediction of the RQKE method systematically improves 
with the expansion order and can fully reproduce the exact quantum dynamics calculated 
from the hierarchy equation. With specihc parameters considered in this paper, the time- 
integrated RQKE rate at the sixth order can be almost identical to the exact result for both 
unbiased and biased system, with both weak and strong site-site coupling strengths. More 
importantly, the detailed time evolution can be exactly predicted as well, as higher-order 
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correction terms are gradually included. The temperature dependence of the equilibrium 
population is also verified, as the classical Boltzmann distribution of the second-order FGR 
prediction is improved toward the exact quantum Boltzmann distribution. The convergence 
order generally increases with the increase of the site-site coupling strength, the decrease of 
the reorganization energy and the decrease of temperature. 

The numerical calculations of this paper are focused on the harmonic bath with a quantum 
Debye spectral density. The formal expression of the QKE rate kernel in Eq. (J5J) is however 
invariant of the bath structure, whether Gaussian or non-Gaussian, so that the RQKE 
method can be applied to a general bath, combined with other numerical methods. The 
mathematical strategy of applying the continued fraction form is not limited to the two- 
site system, and its application to more complicated systems will be demonstrated in our 
forthcoming papers. The RQKE method provides a systematically converged approach of 
quantum dynamics, and its continued fraction form can inspire possibilities of other higher- 
order resummation techniques, such as the extension of the Landau-Zener approximation 
and modifications originally for the lowest order correction. 
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Appendix A: Fourth- and Sixth-Order Quantum Rate Kernels in the two-Site Sys¬ 
tem 


In this appendix, we summarize the expressions of the fourth- and sixth-order QKE rate 
kernels in the two-site system with a 5-spatial correlation. Notice that such a two-site 
system coupled with the harmonic bath is equivalent to the standard spin-boson model with 
a doubled reorganization energy. The fourth-order QKE rate kernel for a general multi-site 
system is derived in Ref. 16], and we simplify this expression with the consideration of the 
two-site system. The forward transfer rate kernel from the donor site 1 to the acceptor site 
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2 is written explicitly as 

^2i ) ( T 2,r3,r 4 )/2|J| 4 

= Re{e iil2T * ~ 2S * ( e 2F * - 1) + e iil2T *~ 2G * ( e 2F4+ - 1) 

+ e lil2T *~ 2G * (e~ 2F * - 1) + e l£l2r t~ 2g ^ (e _2i?4 - 1)}, (Al) 

with rf = r 2 ± r 4 , Gf = g(r 2 ) + s(±t 4 ), and = g(±r 3 ) - g{r 2 + r 3 ) - ^(±(r 3 + r 4 )) + 
g{i 2 + t 3 + r 4 ). 

The sixth-order quantum rate kernel after expanding each term is given by 

/C (6) (r 2 ,--- ,t 6 ) 

= Tr B {F(r 6 )Wp(r 5 )F(r 4 )W P (r 3 )F(r 2 )P e ^} 
~^ 2 \t^ 2 \t a )^ 2 \t 2 ) 

+/C (2) (r 6 )/C (2) (r 4 , r 3 , r 2 ) 

+/C (4) (r 6 ,r 5 ,r 4 )/C (2) (r 2 ). (A2) 

For conciseness, we only present one off-diagonal element of 

y = Tr iJ {K(T 6 )Mp(T 5 )K(T 1 )Mp(T 3 )'R(T 2 )pl;l} 

-/C< 2 )(T 6 )Kl 2 >(T4)V 2 >(T 2 ), (A3) 

and all the other terms can be found from the second- and fourth-order QKE rate kernels. For 
the quantum transport process from the donor site 1 to the acceptor site 2, the corresponding 
term y 2 i is explicitly given by 

3V(-2|J| 6 ) 

= Re jg 4£12r 6 ++ ~ 2 ^ + (g' F 6A + - 1) + e iei2Tfi + ~ 2G 6 + ( e - F <fr> _ i) 

_|_ e iei2T+~-2S++ ( 6 -FqX~ — 1) + e i£ ^ T 6~~ 2g 6 + ( e ~ F 6D + — i) 

-|_ e^ 12T 6 ~ 2g 6 (g ~ F 6B + — 1) + e* ?12T6 ~ 2g 6 fg ~ F 6C — 

_|_ g*ei2'T 6 l_+ —20g —Fqb — 1) + g*^ 12T 6 +— 2g & (g ~ F 6C + — 1^) 

-f e *£i 2 Tg — 26 g ((>~ f 6a + — 1) + e iilaT *~ ~ 2S *~ (e~ F M°~ — 1) 

+ g^ 12 T 6 ++ “ 2 ^ (g ~ F 6A — 1) + g ®^ 127 6 + ~ 2g 6 (g ~ F 6D + — 4) 

_|_ g*ei2r 6 K+ -2CJ 6 + ^g-i ? 6 ^ + _ _|_ g* e l2"r 6 + ~2G 6 + (q~ F 6C — 1) 

+e i £ -i2r 6 +--2e 6 -+ (g -F-j- _ 4 ) + g«T 2 r a ---20 6 -+ -F~ 6 + _ ^ | ( A4 ) 
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Here we introduce the abbreviated notations, r^ =± = r 2 ± 74 ± tq, and = g(r 2 ) + 
g(±r 4 ) + g(±Te), where the left and right ± superscript symbols are associated with r 4 and 
r 6 , respectively. Additional abbreviated notations, r tJ = r* + Tj, = t* + Tj At*,, • • • 
(i,j, k = 2, 3, • • • , 6 ), are introduced to express the functions of F 6 as 

F tjt = 2 9 (t 3 ) ± 2^(±r 5 ) - 2^(r 23 ) - 2^(r 34 ) 

= F25 '(t 45 ) =f 2g(±r 56 ) + 2g(r 234 ) 

±2fi'(r 345 ) ± 25f(r 456 ) =f 25f(r 2345 ) 

= F2(?(r 34 5 6 ) ± 2(yf(r 23456 ), (A5a) 

^6S ± = 2 9(-T 3 ) ± 2^(±t 5 ) - 2g(r 23 ) - 25f(-r 34 ) 

=F 25 f(-r 45 ) =F 25 f(±r 56 ) + 2£f(r 234 ) 

=t2g(—r 34 5) ± 2(yf(—r 4 56) =F 2g(r 234 5) 

=F2fi'(-r 3456 ) ± 25f(r 23456 ), (A5b) 

^6C ± = 2^(t 3 ) =f 2^(±t 5 ) + 2g(r 23 ) + 25-(r 34 ) 

±25f(-r 45 ) ± 2(yf(±r 56 ) - 25f(r 234 ) 

=t2(yf(r 34 5) =F 2g(—r 456 ) =f 2g(r 234 5) 

= F 25 '(t 3456 ) ± 2(yf(r 23456 ), (A5c) 

F oD^ = -2fi'(-r 3 ) =F 2 #(±t 5 ) + 2g(r 23 ) + 25f(-r 34 ) 

±25f(r 45 ) ± 25f(±r 56 ) - 2g{r 23 A) 

±2g(— t 34 s) =F 2(?(t 4 56) =F 2(yf(r 234 5) 

-F2(7(— r 34 56) ± 2(yf(r 234 56), (A5d) 

where the left ± superscript symbol is associated with operations between g functions, and 
the right ± superscript symbol is associated with the sign of the time variable inside g 
functions. 
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